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Abstract 

' We introduce a model for a sandpile, with N sites, critical height N and each site 

connected to every other site. It is thus a mean-field model in the spin-glass sense. We 



find an exact solution for the steady state probability distribution of avalanche sizes, and 
discuss its asymptotics for large N. 
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The publication of a series of articles by P. Bak and collaborators |2| generated a growing 

o 

interest in the study of certain class of cellular automata models now commonly known as 
"sandpiles" because of a crude analogy between the dynamical rules and the way sand topples 
when building a sandpile. The basic motive was not to model the way real sandpiles behave, 
although it did broaden the audience for granular flow (see e.g. || [|, [|), but to study a 
phenomenon they coined "self-organized criticality" in which a system reaches a critical state, 
i.e. a state with no intrinsic length or time scales, without the tuning of an external parameter. 
This generated a big influx of theoretical, computational and even experimental studies of self- 
organized criticality and of sandpiles in general. Unfortunately, even though the number of 
models (i.e. sets of dynamical rules) exhibiting self-organized criticality or falling into the more 
general category of sandpiles has increased dramatically, very few exact results are known at 
the present (see e.g. || |7|, ||, [| [K|). In this letter we present a model in which many quantities 

can be calculated exactly in a rigorous way. 
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A sandpile model is basically a set of dynamical rules describing the way that grains of sand 
are added to a system, the conditions under which those grains can be redistributed inside the 
system, and the way they are removed from the system. Here we consider a system of N sites 
and define h(i) as the (integer) height of the column of sand at site i, i € {1 . . . N}. We drop a 
grain of sand on a site i chosen at random, thereby increasing its height by one: h(i) — > h(i) + l. 
If this new height exceeds the maximum stable value h c then that column topples and gives 1 
grain of sand to each of the N — 1 other sites while one grains drops out of the system. (We 
take h c > N so that h(i) > 0; in fact we are primarily interested in h c = N.) We then examine 
the system to see if any site has a column exceeding h c in which case we topple that column 
also. We keep toppling until all the sites are stable (this characterizes an avalanche). We then 
repeat the procedure of adding a grain at a randomly chosen site. 

This model falls into the category of abelian sandpiles since we always obtain the same 
stable configuration from an unstable one irrespective of the order in which we executed the 
topplings. Dhar used this fact to obtain several properties of those models which we use 
extensively in our analysis. Furthermore, since all sites are connected to all other sites, we can 
say that we have a mean field theory of abelian sandpiles. It is this combination that permits 
the model to be solved. It should be pointed out that our model is quite different from previous 
attempts to study mean-field sandpiles [^, . 

It is shown in p| that not all configurations are allowed in the asymptotic regime, but 
that the allowed (recurrent) configurations all have equal weight. In general, the number of 
recurrent configuration is given by the determinant of the toppling matrix A which is an N x N 
matrix where An = h c and Aj,- = — 1 for each site j connected to i, so that row i of A represents 
the amount of sand lost by every site when site i topples. Thus in the model considered here 
this matrix takes the form 



where 5ij is the Kronecker delta. 

Because of the highly symmetric nature of A in our model it is straightforward to calculate 
the determinant and determine that the number of recurrent configurations is 



A« 



1 + (h e + !)<%, 




Z(N, h c ) = det A = (h c + 1 - N)(h c + 1) 



N-l 
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For the case h c = N we have Z(N,N) = (N + 1) . Incidentally this is the number of 
spanning trees on a fully connected graph with N + 1 sites. Such a result was expected from 
the general observation that det A is precisely the cofactor of the matrix tree (see e.g. fbll ) for 
graphs defined on a superset of the sandpile lattice (there is an additional site, the ground, 
connected to every site). In fact there is a one-to-one correspondence between the recurrent 
configurations and the spanning trees, in the sense that a configuration is recurrent if and 
only if it passes the "burning algorithm" test described in || Actually this holds only for 
symmetric abelian sandpiles, i.e. models where A is symmetric; for the asymmetric case see 




Clearly the recurrent configurations fall into equivalent classes, where two configurations 
are equivalent if one is a permutation of the other. Because of the symmetry of the model 
we need only examine one member of each equivalence class; we find it convenient to restrict 
ourselves to the case h(i) < h[i + 1), i = 1 . . . N — 1, where the amount of sand increases as 
the site label increases. 

Our analysis will focus on a particular subset of the recurrent configurations, namely the 
minimal configurations. A given configuration is minimal if there is no recurrent configuration 
that can be obtained by removing sand from the given configuration. Some brief computations 
show (one can use the so-called "burning algorithm" of §) that the configuration h(i) = i is 
a minimal configuration. All minimal configurations fall into this equivalence class and thus 
there are exactly AH minimal configurations, the permutations of h(i) = i. If we now consider 
equivalence classes for general configurations, we see that the restriction that a recurrent 
configuration be greater than or equal to the minimal configuration yields 

h(i) > i, h(i) > h(i - 1) (3) 

(set h(0) = 0). 

Now consider how avalanches of various sizes come about in our model; we limit the 
discussion to the case h c = N. Size here means either the number of grains that fall out 
of the system or the number of sites that topple; in the model considered here there is no 
ambiguity because they are equal. An avalanche begins when a grain of sand is dropped on a 
site of height N. It will continue if the second-highest site is also at height N. It will be at 
least of size three if the third-highest site is at height N — 1, and it will be at least of size four if 
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the fourth-highest site is at height N — 2, etc. Therefore the size will be (using the equivalence 
class notation) 

Naval = N - min{j : h(i) > i for all j < i < N}. (4) 

This is the only non-zero size avalanche this configuration can produce, and it is determined 
simply by "how long" the configuration stays away from a minimal one. 

If we have such a configuration, what is the probability of starting an avalanche? It is simply 
the fraction of sites at the maximum value, \{i : h(i) = N}\/N. Therefore, to determine the 
avalanche size probability distribution, we merely have to count the configurations with fixed 
iV ava i and with a fixed number of sites at height N. Luckily, the number of such configurations 
can be simply represented as partition functions for smaller systems with different critical 
heights — we obtain (for 1 < k < N): 

P(avalanche of size k) = 

Z(N, N)-'Z(N -k,N-k)j2 ( k _* N _ k ) ~ * * " 2 )- (5) 

A brief explanation of the source of the terms in (|5|) is in order. First of all, since site 
N — k must be at height N — k to get an avalanche of the right size, the remaining sites 
below N — k yield a factor Z(N — k,N — k). Next j is the number of sites at height N. The 
sites between N — k + 1 and N — j must be greater than minimal yet less than N; this is 
equivalent to Z(k — j, k — 2). The combinatoric factor gives us the number of ways we can 
choose configurations with fixed equivalence classes for the subsystems. Finally we have the 
probability of toppling j/N and the normalization Z(N, N). 

Simpler expressions hold for k = 0, 1: 

P(avalanche of size 1) = Z(N — 1,N— 1)/Z(N, N); (6) 

P(no avalanche) = ^ I . -—^-Z(N -j,N- 1) /Z(N, N). (7) 
j=l V J J N 

The above sums (§-0) can be performed exactly. The result is 

Pfavalanche of size k) - (N - k + l) N -"-^-\N 

P(avalanche of size k) - {N + 1)N -i {k _ 1)[{N _ k)[ ( 8 ) 



for 1 < k < N, and 



N-l 

P(no avalanche) = j-j—j- (9) 



The probability distribution (n) has the unusual feature that (for k 7^ 0) it is symmetric about 
fc = (JV + l)/2. 

We are interested in the large N behavior of this result. Taking k > fixed and N — > oo, 
(|) yields 

P(avalanche of size k) ~ — jyp (10) 
and if k is also large (1 <C k <C iV) then this reduces to 

k -3/2 

i-*(avalanche of size k) ~ , — . (11) 



Thus we reproduce the exponent —3/2 previously derived on trees and via numerical and 



nonrigorous arguments fig, |l7j. 



Another quantity of interest is the single-site probability distribution for the heights, i.e. 
the probability that a given site contains H grains of sand in the stationary state. First look 
at configurations where exactly K sites have height H. (In the equivalence class notation 
these sites must be consecutive.) Remove these sites and consider the subsystem consisting 
of the remaining N — K sites. The recurrent configurations of this subsystem are in one-to- 
one correspondence with those of a system of size N — K where the recurrent configurations 
are determined by, in addition to the usual occupancy restrictions (see (||)), the conditions 
j < h(J) < N - 1 if j < H - K and j + K - 1 < h(j) < N - 1 if j > H - K. 

The direct evaluation of the number of allowed configurations satisfying the above criteria 
is possible but somewhat complicated; appropriately summing over K would yield the desired 
probabilities. However, it turns out to be simpler instead to compare the number of allowed 
configurations with K sites having H grains with the number of allowed configurations with K 
sites having H — 1 grains, e.g. by subtracting the latter from the prior. Which configurations 
are left? 

In both cases j < h(j) < N-l for j < H-K and j+K-1 < h(J) < N-l if j > H-K. The 
only difference occurs at site H — K which can be occupied by between H — 1 and N — 1 grains 
in the first case, and between H — K and iV — 1 grains in the second. If K = 1 these conditions 
are identical and no allowed configurations are left. For K > 2, however, site H — K is now 
restricted so that h(H — K) < H — 2. Since we need h(j + 1) > h(j) (for the configuration to 
be a recurrent one), this implies that the first H — K sites can have a maximum of H — 2 grains 



each and thus contribute a factor Z{H — K, H — 2) to the partition function, excluding (global) 
symmetry factors. On the other hand, (|3|) imposes no additional conditions on the last N — H 
sites, and it is easily seen that they contribute a factor Z(N — H, N — H). Defining P{H) = the 
probability that a given site will have H grains, putting in the appropriate symmetry factors 
and summing over the index K we get 

H I nj \ jv- 

P{H) _ P{H _ i) = Z (N, N)- l Z{N-H,N-H)Y,[ KH _ KN _ H \^Z{H-K,H-2), 

(12) 

which is exactly (||)! Since -P(O) = we can collapse the telescoping sum and obtain 

H 

P(site has H grains) = P(avalanche of size k). (13) 

k=l 

Because of the symmetry of the avalanche distribution, we have here the symmetry 
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P(site has H grains) + P(site has N — H grains) = — — - (14) 

for < H < N. 

The asymptotics of (|l3|) are relatively easy to compute using our previous results. For 
N ^> 1 we have 

P(site has H grains) ~ — > — (15) 

V ^ iV^(fc-l)! 17 

which can be used directly to compute the behavior for small H. For N 3> H 3> 1 we have 



P(site has H grains) ~ Jf^l ~ \j ^ J ■ ( 16 ) 

We can also examine the total amount of sand in the system. The distribution of masses 
(total number of grains of sand) of recurrent configurations is most easily computed for mass 
M near the maximum value iV 2 . For example, the number of recurrent configurations with 
mass A^ 2 — 2 is simply N(N + 1)/2. This is easily generalized to arbitrary mass (using inclusion- 
exclusion arguments) and we find that the number of recurrent configurations with total mass 
N 2 -K, Z N 2_ K (N,N), is 



N — 1 



where 

l(N, K)=\n~ (N(N -S) + 2K+ 17/4) 1 / 2 - 3/2] (18) 

simply indicates "how far" one must go in the inclusion-exclusion; \A\ is the least integer > A. 
Then the probability of a configuration having mass M is Zm{N, A^)/(A^+l) iV_1 . Unfortunately 



we are only able to find the asymptotic behavior of (17) for the tails of the distribution, i.e. 
M near N{N + l)/2 or M near TV 2 , which is not where most of the (recurrent) configurations 
reside. Expressions similar to that which were used to find the avalanche distribution, e.g. (||), 
are available: 

N 

: \ \ 

(19) 



Z M (N, he) = £ Z M- jhc (N ~ j, h c - 1) (^j . 



r- 

and one could consider a generating function approach to eliminate the restrictions on j in 
the sum. Unfortunately the resulting expression will only converge in the regions equivalent to 
very large or very small mass, and thus provide no new information. 

The reformulation of the iV-dimensional model in terms of a one dimensional sandpile- 
like problem (with the required symmetry factor attached to each configuration) greatly 
facilitated the computations of "static" quantities like the avalanche size, single-site height, 
and total mass distributions. In order to study the "dynamical" quantities associated to 
the evolution of an avalanche, another reformulation in term of a one dimensional particle 
model proves to be useful; we will use it to calculate the distribution of avalanche durations. 
This reformulation also provides a much more efficient method for numerical study than a 
straightforward implementation of the sandpile process. 

The duration of an avalanche is the number of sweeps that are required in order to reach 
a stable configuration once a grain of sand is dropped on the system. A sweep consists of two 
steps: we first go through all the sites and mark those with a number of grains exceeding the 
critical value and then topple all those sites simultaneously. The process is repeated until no 
site can topple. Note that this definition of duration is not the only one; others are possible 
because of the abelian nature of the model. 

The particle model consists of iV particles moving on a one dimensional lattice with sites 
labeled from left to right by k £ [1, iV + 1] (where N is the number of site in the original model). 
The correspondence between sandpile configurations and particle configurations is as follows: 
given a sandpile configuration Co place a particle at position k for every site in Cq containing k 



grains. If n/c is the resulting number of particles at site k, the stability condition (|3|) translates 
to < k and njsr > 1; since the critical height is iV we also have n^r+i = 0. 

The boundary conditions are "almost" periodic as will be explained below; we also introduce 
a marker (barrier) between sites N and N + 1 that will play an important role in the dynamics. 

The dynamics of the particle model is as follow: pick a particle at random and move it one 
site to the right (this is equivalent of dropping a grain of a sand on Co). If the barrier was not 
crossed then we have a new stable configuration and can repeat the step. Thus, neglecting the 
boundary, our model is a zero-range process; however, the boundary plays a very important 
role. 

If the chosen particle jumps through the barrier it means that one site in Co is unstable and 
will topple. In the particle model such a toppling will mean that all particles simultaneously 
jump one site to the right except for the one at site N + 1 which should move back to site 1. 
An equivalent formulation which we find more convenient is to move the barrier to the left by 
one unit and relabel the sites accordingly (i.e. the barrier defines the position of site N and 
N+l on the periodic lattice). 

As we moved the barrier, a number of particles (say p 2 ) might have gone through it. In the 
sandpile model this means that the initial toppling caused p 2 more sites to become unstable. To 
duplicate their toppling the barrier should now be moved p2 more units to the left. The process 
continues until no particle crosses the barrier as it jumps. The total number of barrier jumps 
required is the duration of the avalanche. The probability that an avalanche has duration T is 
equivalent to the fraction of all particle configurations that will make the barrier jump T times 
if a particle at site N is picked. Of course T = corresponds to the no avalanche case so that 
the probability that the avalanche duration is zero is (N — 1)/(N + 1) (see (Q)). 

If there are T barrier jumps, let Pi,p%, ■ ■ - Pt be the number of particles crossing the 
boundary at each jump and let P = p\ + p2 + ■ ■ ■ + Pt] note that pi > for i < T and 
Pt = 0; we also have p\ > 1 for T > 1. Provided T > we have 

P(durationT) = Z(N,N)- 1 ]T ^W N (l; Pu p 2 , . . . ,p T ) (20) 

pi,...,PT restricted 

where Wjsr(l;pi,p2, ■ ■ ■ ,Pt) is the number of configurations that produce the specified barrier 
jumps. We can compute this in an iterative fashion, in terms of WM(s',Pi,Pi+i, ■ ■ ■ ,Pt), the 
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number of (restricted) configurations on the subsystem consisting of sites 1 . . . M; s is size of 
the current jump, and pi,Pi+±, ■ ■ ■ ,pr are the remaining jump sizes. We have 

[N\ 

W n (1;pi,P2,-..,Pt) = \Jw(p 1 ,l)W N - 1 (p 1 -l;p2,P3,---,PT), (21) 

[n _ Pl \ 

W N - 1 (p 1 -1;p2,P3,...,pt) = \w(p2,pi-l)W N - pi (p 2 ;p3,P4,...,PT), (22) 

V P2 J 

and 

W M (pi-i;Pi,Pi+i,---,Pr) = (^^w{Pi,Pi-i)W M - Pi _ 1 (Pi;Pi+i,Pi+2,---,PT) (23) 

for 2 < % < T — 1. The weight w(p, s) is easily computed since the criteria for allowed 
configurations are automatically satisfied within each block. Since each of the p particles can 
independently be at any of the s sites, 

w(p, s) = s p . (24) 

Finally, the last step of the iteration is simply 

(N - p + VT A 

WN-p+rr^ipr-KPr) = U;(p T ,p T _i)Z(iV - P, N - P) = Z(N - P, N - P). 

\ PT J 

(25) 

As an example, for T = 4 this would lead to 

PM ,• 1! ^2 N-g-1 N-n-n Z (N-P,N-P) n ( N \ 
F(d„ra. 1 cn4)=S £ S 2( „, iV ( P1 KPspJ ^ " " ^ ' P6) 

The general case can actually be written in a reasonably compact form. We can rearrange 
the terms and obtain the following (for T > 0): 

TV 



5' 

I St! 



(27) 

where we have the restriction si = 1, Sj > 0. We identify Sj with except that because of 
the need to initiate the avalanche s\ = 1 and s 2 = pi — 1- Of course is nothing more than 
the size of the i th jump. 
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